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1  INTRODUCTION 


The  use  of  computer  simulations  as  an  optimal  design  tool  which  lessens  the  costs  in  time 
and  effort  in  experimental  and  physical  testing  has  become  prevalent  in  aerospace  design 
( e.g .,  the  design  of  the  Boeing  777)  and  development  of  automotive  mass  production 
items.  It  is  also  commonly  found  in  the  design  of  production  and  assembly  lines  where 
repetitive  motion  related  to  parts  assembly  is  critical.  In  this  presentation  we  discuss  a 
non-aerospace,  noil-automotive  design  use  of  computer  simulations. 

A  team  at  N.C.  State  University  composed  of  material  scientists,  physicists,  and 
applied  mathematicians,  have  used  computer  simulations  as  a  fundamental  design  tool 
in  developing  a  new  prototype  high  pressure  organometallic  chemical  vapor  deposition 
(HPOMCVD)  reactor  for  use  in  thin  him  crystal  growth.  The  advantages  of  such  a 
reactor  lie  in  improved  advanced  technology  products  involving  exotic  materials  with 
high  thermal  decomposition  pressures  and  increased  control  over  local  stoichiometry  and 
defect  formation. 

While  we  focus  here  on  reactor  design  geometry,  innovative  use  of  computer  simulation 
as  a  design  tool  has  also  been  made  with  regard  to  sensing  -  development  of  a  polarized 
reflectance  spectroscopy  (PRS)  real  time  monitoring  system  (e.g..  [2,  3])  -  as  well  as 
innovative  nonlinear  feedback  control  design.  Considered  with  the  work  we  describe 
below,  this  effort  offers  a  strong  endorsement  of  such  multidisciplinary,  computationally 
based  modeling  teams  in  the  design  of  new  products  in  areas  of  emerging  technologies 
where  heretofore  extensive  and  costly  experimental  design  was  the  central  paradigm. 


1.1  High  Pressure  OMCVD 

A  simplified  view  of  the  organometallic  chemical  vapor  deposition  (OMCVD)  process 
is  depicted  in  Fig.  1  for  the  case  of  gallium  phosphide  him  growth.  For  OMCVD  at 
elevated  pressures  the  source  material  (trimethylgallium  and  tertiarybutylphosphine)  is 
transported  to  the  hot  substrate  by  a  carrier  gas  (nitrogen).  When  exposed  to  elevated 
temperatures,  the  source  gas  pyrolytically  decomposes.  Molecular  fragments  and  atoms 
that  diffuse  to  the  hot  substrate  stick  to  the  surface  where  they  diffuse,  undergo  further 
chemical  reactions,  and  are  incorporated  into  the  growing  layer  (GaP).  Operations  at 
elevated  pressures  would  broaden  the  applicability  of  OMCVD  to  include  materials  that 
exhibit  high  decomposition  pressures  and  would  allow  increase  in  the  growth  tempera¬ 
tures  for  materials  presently  used. 

For  heteroepitaxy,  source  gases  may  be  transported  in  alternating  pulses  as  alter- 
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Figure  1:  Schematic  depiction  of  organometallic  chemical  vapor  deposition  at  elevated  pressure. 


nating  layers  are  grown  on  the  substrate.  Optimal  design  of  a  high  pressure  OMCVD 
reactor  requires  uniform  and  efficient  transport  of  source  material  to  the  hot  substrate, 
without  excessive  retention  times  of  the  source  material,  in  order  to  achieve  uniform  him 
thicknesses  with  abrupt  transitions  from  one  layer  to  the  next.  Because  of  these  require¬ 
ments,  special  care  must  be  taken  to  avoid  formation  of  recirculation  cell  regions,  which 
are  more  likely  to  occur  at  elevated  pressures  and  temperatures. 


2  OPTIMAL  DESIGN  SIMULATIONS 


The  optimal  design  process  utilized  computer  simulation  as  a  design  tool.  As  a  design 
tool  the  simulations  allowed  rapid  exploration  of  design  choices  and  helped  to  identify 
key  design  issues.  A  simulation  that  includes  every  possible  process  and  reaction  was  not 
feasible.  Rather,  the  aim  of  the  simulations  was  to  provide  insight  into  the  important 
reactor  design  features  and  processes  that  critically  effect  HPOMCVD  growth. 

In  our  investigations  we  employed  both  steady-state  and  time-dependent  simula¬ 
tions  assuming  compressible  flow,  temperature-dependent  density  variations,  and  buoy- 
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ant  forces.  The  governing  equations  were  those  for  mass 

+  V  •  (,m)  =  0, 

momentum 

p(  ff  +  u-Vu)  =  -VP  +  V-T-/9g 

T  =  -|//(V  ■  u)I  +  //(Vu  +  VuT), 

energy 

^ll  +  u-Vhj  =  V  •  (A’VT), 
and  the  A4h  gas  phase  species 

P  “1“  u  '  =  ^  '  {pDk^Yk)  +  Rky 

where  g  is  the  gravitational  acceleration,  u  is  the  velocity,  T  is  the  temperature,  Y);  is  the 
mass  fraction  of  the  At  li  chemical  species,  R;.  represents  the  chemical  reaction  rate,  and 
P  is  the  pressure.  The  density  variations  were  modeled  according  to  the  relationship, 
p  =  po [1  —  13{T  —  To)],  where  T0  is  a  reference  temperature,  p0  is  a  reference  density 
calculated  from  the  ideal  gas  law  at  the  reference  temperature  and  reactor  pressure,  and 
ii  is  the  volume  coefficient  of  expansion  ( 0  =  1  /T). 

The  gas  parameters,  dynamic  viscosity  (//),  specific  heat  (cp),  and  conductivity  (A’), 
and  the  gas  phase  diffusivity  of  the  species  ( Dk )  have  temperature-dependent  values.  For 
our  simulations,  the  temperature-dependent  values  of  the  gas  parameters  ( p,cp,k )  were 
linearly  interpolated  from  measured  values  taken  from  the  available  literature  [15,  16, 
4],  The  viscosity  and  specific  heat  remain  essentially  unchanged  (<  1%  variation)  for 
pressures  up  to  10  atmospheres  [11],  and  the  conductivity  for  the  case  of  N2  is  similarly 
expected  to  be  independent  of  pressure  for  pressures  up  to  about  10  atmospheres  [7]. 

The  governing  equations  were  discretized  using  the  Galerkin  finite  element  method 
with  weighted  residuals  for  the  degrees  of  freedom.  A  mixed  formulation  with  quadri¬ 
lateral  elements  was  used  with  piecewise  linear  discontinuous  elements  for  pressure,  and 
quadratic  elements  for  the  other  degrees  of  freedom.  Simulations  were  performed  with 
commercially  available  code  (FIDAP,  Fluid  Dynamics  International,  Evanston,  IL)  on  a 
Silicon  Graphics  Power  Indigo  2. 


3  VERTICAL  CYLINDRICAL  DESIGN 


Because  of  their  inherent  strength,  we  initially  focused  on  the  evaluation  of  reactors  of 
cylindrical  geometry,  with  central  (axial)  injection  of  source-laden  carrier  gas  streams 
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Figure  2:  An  early  design  using  cylindrical  geometry. 


that  are  symmetrically  distributed  into  a  radial  flow  channel  and  exhausted  via  an  an¬ 
nular  exhaust  channel.  The  early  evolution  of  the  cylindrical  reactor  design  was  rapid 
as  simulation  results,  evaluation  of  the  available  literature,  and  design  considerations 
were  assimilated.  An  early  example  of  such  a  design  (Fig.  2)  entailed  a  trumpet-shaped, 
central  injection  nozzle  containing  a  set  of  nested  injection  ports  and  a  single  substrate 
centered  on  the  axis-of-symmetry.  Two  potential  problem  areas,  the  formation  of  recir¬ 
culation  cells  between  the  trumpet-shaped  inlet  nozzle  and  the  outer  cylindrical  wall  and 
the  build-up  of  source  material  at  the  center  of  the  wafer,  were  demonstrated  by  CFD 
simulations  of  this  reactor  geometry.  Retention  of  source  material  in  recirculation  cells 
is  quite  likely  to  degrade  the  efficient  switching  necessary  for  achieving  sharp  transitions 
between  heteroepitaxial  layers  by  allowing  source  material  to  diffuse  to  the  substrate 
after  termination  of  the  source  pulse. 

Fig.  3  depicts  a  later  version  of  the  cylindrical  reactor  design  in  which  the  wafer  has 
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been  moved  off  axis.  This  avoids  the  build-up  of  source  material  at  the  wafer  center 
and  allows  multi-wafer  processing.  Fig.  3  also  reveals  that  the  trumpet-shaped  source 
nozzle  has  been  replaced  by  a  cylindrical  nozzle  with  a  water-cooled  inner  core.  In  the 


Figure  3:  Subsequent  design  with  the  wafer  moved  off-axis. 

final  version  of  the  cylindrical  vertical  reactor  design  (Fig.  4)  the  single  susceptor  of 
the  previous  design  has  been  replaced  by  a  heater  assembly  block  that  contains  separate 
susceptors  for  each  wafer  and  clearly  indicates  the  multiwafer  capability  of  this  design. 
Because  simulations  indicated  that  recirculation  cells  may  form  where  the  channel  flow 
joined  the  purge  flow  at  the  outer  annulus,  the  outer  edge  of  the  heater  assembly  block 
was  shaped  to  guide  the  carrier  gas  around  the  90°  bend,  so  that  carrier  and  purge  gas 
flow  directions  were  parallel  when  they  joined. 

Early  on  in  the  design  process  the  simulations  included  both  solid  and  gas  portions 
of  the  model.  As  the  design  evolved  the  simulations  focused  in  on  the  essential  gas  phase 
components  while  solids  were  modeled  by  applying  appropriate  boundary  conditions  (no- 
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slip,  heat  transfer,  .  .  .).  In  further  design  evolution,  the  transfer  of  the  gas  from  the 
central  injection  tube  to  the  radial  flow  was  identified  via  simulations  to  be  another 
critical  design  feature.  Extensive  simulations  with  the  model  shown  in  Fig.  5  led  to  a  set 
of  design  criteria  for  the  vertical  cylindrical  reactor  design  at  elevated  pressures  [13]. 


W 


Purge  Inlet  y 


Outlet  Inlet 


Figure  5:  Model  geometry  showing  cylindrical  symmetry  about  the  z-axis,  with  R=6A  mm,  dw—  25  mm, 
W= 89  mm,  and  rw  —  51  mm.  The  gravitational  acceleration  was  either  oriented  along  the  positive 
direction  of  the  z-axis  (wafer  up)  or  along  the  negative  direction  of  the  z-axis  (wafer  down).  The  dashed 
line  perpendicular  to  the  substrate  indicates  the  position  of  temperature  and  velocity  profiles  used  in 
subsequent  figures. 


3.1  Atmospheric  pressure  results. 

In  order  to  establish  a  baseline  for  comparison  at  elevated  pressures,  we  first  performed 
simulations  at  atmospheric  pressure.  The  conditions  of  these  simulations  were  as  follows. 
The  channel  height  varied  from  0.8  to  12.7  mm.  The  gravitational  acceleration  pointed 
away  from  wafer  surface.  Inlet  and  purge  flows  were  each  12  standard  liters  per  minute 
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(slpm).  The  following  temperature  profile  was  established  on  the  heater  assembly  block 
surface:  1200°K  on  the  wafer  surface  decreasing  linearly  from  the  edge  of  the  wafer 
in  either  direction.  The  inlet  temperature  was  fixed  to  298°K  by  water  cooling  of  the 
inlet  tube.  Heat  fluxes  were  specified  for  the  thermal  boundary  conditions  at  the  inner 
(lower  channel  wall)  and  outer  walls,  q  =  h[T  —  7), jj,  where  Tref  =  350  and  400°K 
for  the  inner  and  outer  wall,  respectively  (Fig.  5).  Fig.  6  shows  (a)  the  temperature 


Figure  6:  Atmospheric  pressure  simulations  with  12  slpm  carrier  and  purge  flows,  showing  (a)  tempera¬ 
ture  contour  plot  for  1200°K  wafer  temperature,  where  contours  represent  steps  of  180°K  ,  Ti=1110°K  , 
and  Tv =390°  K  and  streamlines  contour  plots  for  (b)  1200°K  and  9. 5-mm  channel  height  (A=4.1e-4  and 
B=-1.5e-4),  (c)  a  room  temperature  simulation  and  9. 5-mm  channel  height  (A=1.8e-4  and  B=4.0e-3), 
and  (d)  1200°K  and  3.2-mm  channel  height  (A=4.3e-4  and  B=-3.4e-4). 
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distribution  in  the  channel  and  (b)  the  flow  streamlines  for  simulations  with  the  above 
conditions.  Because  of  the  water  cooling  at  the  inner  wall  the  heating  of  the  vapor  phase 
was  confined  primarily  to  the  channel  region  in  the  immediate  vicinity  of  the  wafer.  The 
calculated  temperature  of  the  lower  channel  surface  opposite  to  the  wafer  was  351°  Iv 
without  accounting  for  radiative  heat  exchange  and  379°K  if  radiative  heat  exchange  was 
included  in  the  simulation. 

The  largest  velocities  (3  m/s)  occurred  at  the  center  of  the  reactor  along  the  axis-of- 
symmetry,  where  the  flow  cross  section  was  the  smallest  (diameter  R  in  Fig.  5).  A  region 
of  high  static  pressure  occurred  at  the  first  90°  bend,  along  the  curved  portion  of  the 
heater  assembly  block,  while  a  region  of  low  static  pressure  occurred  farther  downstream 
along  the  inner  wall  of  the  channel.  This  low  pressure  region  coincided  with  a  recirculation 
flow  region  (Fig.  6b).  A  second  recirculation  region  occurred  farther  downstream,  along 
the  heater  assembly  block  surfaces  in  the  vicinity  of  the  wafer.  The  presence  of  the 
recirculation  regions  constricted  the  flow  in  the  channel  and  increased  the  peak  flow 
velocity  locally. 

The  results  of  Fig.  6  illustrate  that  the  recirculation  regions  at  1200°K  (Fig.  6b)  also 
occurred  at  room  temperature  (Fig.  6c),  that  is,  were  not  thermally  induced,  and  were 
reduced  by  decreasing  the  channel  height  (Fig.  6d).  The  recirculation  regions  appeared 
to  be  caused  by  the  ninety-degree  bend  and  expanding  cross  section  at  the  bend.  This  is 
analogous  to  the  situation  encountered  in  a  diffuser,  where  the  gas  velocity  decreases  as 
the  duct  expands  producing  an  increase  in  static  pressure.  This  adverse  pressure  gradient 
may  cause  separation  of  the  boundary  layers  from  the  diffuser  walls  and  accompanying 
recirculation  regions  [8]. 

The  extent  of  this  diffuser  stall  increases  as  the  ratio  of  the  entrance  to  exit  cross- 
sectional  areas  increases.  The  effect  of  the  channel  height  on  the  recirculation  regions 
can  be  seen  in  Fig.  7,  where  the  radial  velocity  profile  across  the  channel  is  plotted  for 
channel  heights  of  0.8-12.7  mm.  From  Fig.  7  it  can  be  seen  that  at  constant  flow  rate 
and  increasing  channel  height,  i,  the  recirculation  regions  increasingly  confine  the  bulk 
flow  to  a  small  portion  of  the  channel  area.  At  a  value  of  L  where  the  cross-sectional 
area  after  the  bend  equals  the  cross-sectional  area  before  the  bend  (Fig.  7  curve  b), 
the  flow  profile  takes  on  the  parabolic-type  distribution  expected  for  channel  flow.  Thus 
in  reactor  designs  necessitating  a  bend  in  the  inlet  flow  there  exists  a  criterion  for  the 
optimization  of  flow  relating  the  diameter  of  the  inlet  tube  to  the  channel  height  L. 
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Figure  7:  Plot  of  the  radial  component  of  the  velocity  (ur)  across  the  channel  (z)  at  a  point  after  the 
90°  bend.  The  radial  velocity  is  non-dimensional  with  respect  to  the  peak  radial  velocity  (umax),  while 
the  distance  across  the  channel  is  non-dimensional  with  respect  to  the  channel  height  (L).  Curves  are 
plotted  for  channel  height  and  peak  velocities  of  (a)  0.8  mm  and  6.4  m/s,  (b)  1.6  mm  and  3.2  m/s,  (c) 
3.2  mm  and  2.4  m/s,  (d)  6.4  mm  and  2.4  m/s,  (e)  9.6  mm  and  2.4  m/s,  and  (f)  12.7  mm  and  2.4  m/s. 

3.2  Elevated  pressure  results 


Fig.  8  represents  the  results  of  simulations  for  (a)  9.5-mm  and  (b)  3.2-mm  channel 
heights,  respectively,  at  10  atmosphere  pressure.  All  other  conditions  remained  the  same, 
including  the  mass  flow  rate  of  12  slpm  into  the  reactor  for  both  the  carrier  and  purge 
gases.  All  inlet  velocities  were  scaled  accordingly.  As  a  result,  the  Reynolds  number 
remained  unchanged,  because  the  increase  in  density  was  offset  by  the  decrease  in  velocity. 
In  Fig.  8c  we  also  show  the  result  of  a  simulation  for  a  9.5-mm  channel  height,  but 
inverted  geometry,  that  is,  with  gravity  vector  g  pointing  toward  the  wafer  surface.  In 
both  orientations  the  recirculation  region  filled  the  entire  channel,  but  in  the  configuration 
of  Fig.  8c  the  recirculation  region  moved  further  downstream,  that  is,  directly  above  the 
wafer,  and  the  temperature  gradient  from  the  wafer  to  the  opposite  wall  became  non- 
uniform.  For  the  3.2-mm  channel  height,  there  were  no  recirculation  regions  in  the 
vicinity  of  the  substrate  and  the  temperature  gradient  remained  uniform  from  the  wafer 
surface  to  the  surface  of  the  opposing  channel  wall,  for  both  the  normal  and  the  inverted 
(not  shown)  configurations. 
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Figure  8:  Streamline  (left)  and  temperature  (right)  contour  plots  at  10  atmosphere  with  12  slprn 
carrier  and  purge  flows  and  wafer  temperature  of  1200°K  for  (a)  9.5-mm  channel  height  (A=-1.9e-5 
and  B=3.7e-5),  (b)  3. 2-mrn  channel  height  (A=1.9e-4  and  B=5.0e-5),  and  (c)  9.5-mm  channel  height  in 
an  inverted  orientation  (wafer  down)  (A=l.le-4  and  B=-1.8e-4),  where  temperature  contours  represent 
steps  of  180°K  T1=1110°K  ,  and  T2=390°K  . 

3.3  Radiation  and  Opposite  Wall  Deposition  Effects 


Constraining  the  radial  flow  in  the  small  channel  will  act  to  suppress  buoyancy-induced 
recirculation,  which  is  enhanced  at  higher  pressures.  The  channel  heights  here  (0.8- 
12.7  mm)  are  lower  than  those  reported  on  previously  [6,  9,  10,  17]  and,  for  such  small 
channel  heights,  heating  of  the  opposite  inner  wall  becomes  a  concern.  We  investigated 
the  effect  of  radiative  heat  exchange  between  the  hot  graphite  substrate  and  the  water- 
cooled  graphite  inner  wall.  For  the  case  of  the  3.2-mm  channel  height  (1200°K  wafer 
temperature,  12  slpm  flow  rate,  and  atmospheric  pressure)  the  simulated  inner  wall  tem¬ 
perature  opposite  the  hot  substrate  reached  a  temperature  of  379°K  ,  an  increase  of  28°K 
from  a  similar  simulation  without  radiative  heat  exchange.  Even  with  radiative  heating, 
the  inner  wall  temperature  would  still  be  lower  than,  e.g.,  the  pyrolysis  temperature  of 
trimethylgallium  (723  °K  ),  a  widely  used  gallium  source  compound.  A  comparison  of  the 
channel  temperature  profile  with  and  without  radiation  shows  little  difference  between 
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the  two  cases,  except  in  the  vicinity  (0.5  mm)  of  the  inner  wall  (Fig.  9). 

Because  of  the  small  channel  heights  necessary  for  super- atmospheric  OMCVD,  dif¬ 
fusion  of  products  of  the  pyrolysis  of  source  vapors,  such  as  trimethylgallium,  from  high 
temperature  regions  to  the  inner  wall  is  a  concern.  Molecular  fragments  that  diffuse  to 
the  cold  inner  wall  would  deposit  there.  The  formation  of  deposits  can  influence  the  ra¬ 
diative  properties  of  the  inner  wall,  though  this  effect  may  act  to  either  reduce  or  increase 
the  wall  temperature,  depending  upon  the  film  thickness  and  the  wall  material  [5].  The 
importance  of  opposite  wall  deposition  can  be  estimated  by  comparing  the  diffusion  time 
of  source  material  diffusing  to  the  inner  wall  with  the  transit  time  through  the  portion  of 
the  reactor  containing  the  hot  substrate.  These  characteristic  times  will  be  affected  by 
the  dimensions  of  the  channel  height  L  and  the  diameter  of  the  hot  substrate  dw  as  well  as 
the  temperature  profile  in  the  channel  and  the  carrier  gas  flow  rate.  A  diffusion  length  can 
be  defined  as  the  perpendicular  distance  from  the  inner  wall  to  the  location  in  the  chan¬ 
nel  at  which  the  temperature  reaches  the  pyrolysis  temperature  of  TMG  (723°K  )  [14], 
For  the  simulation  illustrated  in  Fig.  9,  this  corresponds  to  a  diffusion  length  of  approx¬ 
imately  1.2  mm.  As  seen  in  Fig.  9  the  temperature  distribution  across  the  channel  is 
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Figure  9:  Plot  of  the  temperature  profile  across  the  channel,  position  as  indicated  by  the  dashed  line 
in  Fig.  5,  for  simulations  with  (x)  and  without  (o)  radiative  heat  exchange  between  the  substrate  and 
the  inner  wall. 

roughly  linear.  Assuming  a  linear  temperature  profile,  the  characteristic  diffusion  time 
for  source  material  to  reach  the  inner  wall  is  G;//  =  ( Tp  —  TW)2L2 /(D(TS  —  Tw )2),  where 
Tp,  7',.,,,  and  Ts  are  the  pyrolysis,  inner  wall,  and  substrate  temperatures,  respectively, 


.5  0  0.5  1  1.5  2  2.5  3  3.5 

Distance  across  the  channel  [mm] 
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and  D  is  the  diffusion  coefficient.  The  characteristic  transit  time  of  material  past  the 
hot  substrate  can  be  estimated  as  transit  =  dw/u{8 ),  where  dw  is  the  substrate  diameter 
and  u(5)  is  the  radial  velocity  one  diffusion  length  from  the  inner  wall.  For  small  enough 
channel  heights,  the  radial  velocity  profile  can  be  assumed  to  be  parabolic  (Fig.  7).  To 
avoid  deposition  of  source  material  on  the  inner  wall  opposite  the  substrate,  the  transit 
time  of  source  material  past  the  region  of  the  hot  substrate  should  be  much  less  than 
the  diffusion  time  of  the  source  material  to  the  inner  wall.  In  Fig.  10  the  ratio  of  the 


Figure  10:  Ratio  of  transit  time  of  source  material  through  the  region  above  the  substrate  to  the 
diffusion  time  of  source  material  to  the  inner  wall  as  a  function  of  the  dimensionless  channel  height 
(channel  height/wafer  diameter),  assuming  the  same  flow  rate  for  all  channel  heights.  The  data  plotted 
here  is  for  the  case  of  TMG  (723°K  pyrolysis  temperature)  in  Nj.  The  deposition  of  source  material  on 
the  inner  wall  will  be  minimized  for  small  ratios  of  transit  to  diffusion  times.  The  dashed  line  indicates 
a  ratio  of  one. 

transit  time  to  the  diffusion  time  are  plotted  as  a  function  of  the  dimensionless  channel 
height,  L/dw ,  for  a  fixed  flow  rate.  The  portion  of  the  curve  lying  below  the  dotted  line 
corresponds  to  ratios  of  less  than  one.  From  Fig.  10  it  is  apparent  that  deposition  of 
source  material  on  the  inner  wall  opposite  to  the  substrate  can  be  avoided  by  choosing  a 
suitable  ratio  of  channel  wall  height  to  wafer  diameter,  either  by  increasing  the  channel 
height  or  decreasing  the  wafer  diameter. 

For  flow  in  rectangular  and  cylindrical  cavities,  forced  convection  is  expected  to  dom¬ 
inate  over  buoyancy-induced  reentrant  flow  when  ( Gr/Re 2)  <<  1  [12].  The  Grashof 
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number,  Gr  =  p2 f3gATL3 / p2 ,  is  a  measure  of  the  ratio  of  buoyancy  to  viscous  forces, 
while  the  Reynold’s  number,  Re  —  pULj p,  is  a  measure  of  the  ratio  of  inertial  to  viscous 
forces,  where  U  is  the  characteristic  velocity  and  p  and  p  have  been  previously  defined. 
Here,  the  characteristic  length  L  is  the  channel  height,  and  AT  is  the  temperature  dif¬ 
ference  across  the  channel.  In  the  vertical  cylindrical  reactor  the  gas  flows  radially  out 
between  the  heater  assembly  block  and  the  inner  (water-cooled)  wall  (Fig.  5).  Therefore, 
AT  and  the  characteristic  velocity,  U,  depend  upon  the  radial  position  in  the  channel, 
since  the  temperature  of  the  heater  assembly  block  has  a  specified  radial  dependence  and 
the  velocity  is  decreasing  with  radial  distance  as  the  cross-sectional  area  increases.  In 
addition,  the  flow  velocity  will  also  depend  on  the  temperature  and  will  increase  with 
increasing  temperature  (smaller  density).  Therefore  we  sought  to  determine  the  charac¬ 
teristic  ratio,  Gr /Re2 ,  that  is  indicative  of  a  flow  regime  dominated  by  inertial  forces.  For 
this  vertical  cylindrical  reactor  design  with  constrained  radial  flow,  simulations  showed 
that  Gr/Re2  <  10  was  indicative  of  the  absence  of  recirculation  regions.  For  a  given 
temperature  difference  between  the  hot  wafer  and  the  cold  opposite  wall,  AT,  and  offset 
rw  of  the  wafer  from  the  vertical  axis  of  the  inlet  tube,  the  minimum  flow  rate,  Q.  thus 
is  determined  by  the  choice  of  L.  Also,  there  exists  a  constraint  on  the  radius,  R ,  of  the 
inlet  tube  imposed  by  the  requirement  of  Re  <  1000  for  noil-turbulent  flow.  Combined 
with  the  above  discussed  optimization  of  flow  for  matching  cross-sectional  areas  before 
and  after  the  inlet  bend,  this  determines  the  inlet  bend  radius  of  curvature,  rc. 

For  example,  with  N2  carrier  and  purge  gas,  1200°K  substrate  temperature,  and  350°K 
inner  wall  temperature  Fig.  11  indicates  graphically  the  required  normalized  entrance 
radius,  defined  to  be  ( R  +  rc)/r„,,  as  a  function  of  the  dimensionless  channel  height, 
L/rw.  Two  curves  are  shown  corresponding  to  reactor  operating  pressures  of  1  and  2 
atmospheres.  The  normalized  entrance  radius  has  a  maximum  value  of  1,  or  the  inlet 
bend  will  extend  above  the  wafer.  The  normalized  entrance  radius  should  lie  on  or  below 
the  line  indicated  for  the  given  operation  pressure. 


4  DESIGN  RULES 


Computer  simulations  of  fluid  and  heat  transport  in  a  vertical  reactor  for  organometallic 
chemical  vapor  deposition  at  elevated  pressure  with  radial  channel  flow  outward  from  a 
central  gas  injection  tube  strongly  suggest  the  following  design  constraints: 

1.  Matching  cross-sectional  areas  before  and  after  the  inlet  bend  is  required  for  an 
optimum  flow  profile  at  the  entrance  of  the  channel.  This  imposes  a  constraint 
on  the  inlet  bend  radius  that  depends  also  on  the  choice  of  channel  height  and 
temperature  drop  across  the  channel  at  the  location  of  the  heated  substrate  wafer. 
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Figure  11:  Plot  of  the  inlet  length  parameter  (radius  of  central  injection  tube  plus  radius  of  curvature 
of  the  inlet  bend)  as  a  function  of  the  channel  height,  for  operating  pressures  of  1  and  2  atmospheres. 
Specified  reactor  conditions  include  of  N2  carrier  gas,  1200°K  substrate  temperature,  and  350°K  inner 
wall  temperature.  Both  lengths  are  non-dimensional  with  respect  to  the  radial  position  of  the  substrate 
on  the  heater  assembly  block.  The  dashed  line  indicates  where  the  inlet  length  parameter  is  equal  to 
the  radial  position  of  the  substrate. 


2.  The  channel  height  is  constrained  by  both  the  need  for  imposing  forced  flow  to 
prevent  buoyancy  related  reentrant  flow  in  the  vicinity  of  the  hot  substrate  and 
preventing  diffuser  stall  at  the  channel  entrance.  The  critical  value  of  channel 
height  at  which  reentrant  flow  ensues  decreases  with  increasing  pressure. 

3.  The  critical  value  of  the  channel  height  imposes  a  constraint  on  the  diameter  of 
wafers  that  can  be  processed  without  deposition  of  fragments  of  the  injected  source 
vapor  molecules  on  the  wall  opposite  to  the  substrate  wafer  to  typically  less  than 
ten  times  the  channel  height. 


The  last  constraint  on  the  diameter  of  the  wafers  can  be  relaxed  by  increasing  the 
flow  rate.  However,  increasing  the  flow  rate  at  high  pressures  may  require  operation  in 
the  turbulent  flow  regime  and  it  would  be  necessary  to  investigate  the  effects  this  may 
have  on  mass  transport  to  the  substrate. 
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5  FLEXIBLE  GEOMETRY  DESIGN 


Based  on  the  vertical  reactor  simulations  we  modified  the  reactor  geometry  to  remove  the 
inlet  bend,  resulting  in  a  horizontal  input  and,  thus,  a  horizontal  flow  reactor.  The  first 
prototype  reactor  built  was  a  rectangular  cross-section  reactor  incorporating  a  horizontal, 
gradual  expansion  duct  between  the  source  gas  lines  and  the  main  reactor  body  (Fig.  12). 
The  high  pressure  CVD  reactor  chamber  is  modular  in  design  allowing  for  modifications. 
Two  reactor  chambers  have  been  built,  the  first  for  atmospheric  pressure  validation  of 
simulations,  permitting  data,  collection  methods  that  are  difficult  to  incorporate  into 
an  industrial  reactor  design,  and  the  second  for  crystal  growth  at  elevated  pressures. 
Operation  at  elevated  pressures  requires  an  external  enclosure  of  the  reactor  chamber, 
prohibiting  such  data  collection  techniques,  but  validation  of  the  flow  inside  a  reactor 
with  similar  thermal  profiles  and  geometry  builds  confidence  in  the  simulations.  The 
simulations  of  the  horizontal  reactor  will  be  used  to  provide  input  for  the  construction 
of  a  reduced  order  model  for  real-time  control.  The  simulation  of  growth  kinetics  and 
comparison  to  experimental  data  will  aid  in  the  identification  of  a  reduced  order  surface 
kinetics  model. 


pressure  confinement  chamber 


substrate  loading 


exhaust 


quartzglass 

reactor 


laser  probe 
beam 


process  gas 


reinjection 


Figure  12:  HPOMCVD  reactor  with  flexible  geometry. 
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The  modeling  equations  to  be  employed  are  the  same  as  those  used  in  the  vertical 
reactor  studies.  Since  there  is  (to-date)  no  active  cooling  of  the  reactor  walls,  radiative 
heat  exchange  between  the  inner  surfaces,  external  radiative  cooling,  and  conduction 
of  heat  in  the  plane  of  the  quartz  walls  become  important  components  of  physically 
meaningful  simulations.  Because  of  the  size  of  the  computer  model  the  three  components 
(expansion  duct,  main  reactor  body,  and  exhaust  duct)  have  been  simulated  separately. 


10  ATMOSPHERE 


Figure  13:  Simulations  of  flow  in  the  expansion  tube  at  10  atmosphere  for  different,  flow  rates.  Interior 
cross-sectional  planes  are  shown  to  enhance  visualization  of  the  flow.  The  simulations  indicate  that 
backflow  is  expected  at  9.0  slpm. 

Simulation  of  the  expansion  duct  predicts  diffuser  stall  at  approximately  10  atmo¬ 
sphere  (Fig.  13).  Further  simulation  predicts  that  backflow  conditions  may  be  avoided 
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by  maintaining  a  constant  cross  section  from  the  tubular  gas  injection  port  to  the  main 
body  of  the  reactor [1],  The  predicted  diffuser  stall  may  afford  another  opportunity  to 
test  the  validity  of  the  simulation  model  once  the  reactor  is  operating.  The  modular 
design  of  the  reactor  permits  the  use  of  a  redesigned  expansion  tube  without  discarding 
other  portions  of  the  reactor,  should  experiments  identify  this  as  a  problem. 

Simulations  in  the  main  body  of  the  reactor  predict  a  longitudinal  vortex  when  the 
substrate  is  heated  (Fig.  14),  but  this  is  not  expected  to  significantly  increase  the  reten¬ 
tion  time  of  the  source  material. 


Figure  14:  Simulation  of  the  main  reactor  body  showing  longitudinal  vortices.  Inverted  arrows  on 
the  left  indicate  the  relative  orientation  of  three  sets  of  massless  particles  at  the  inflow,  while  the  traces 
show  the  path  of  the  particles.  The  arrows  on  the  right  indicate  the  relative  orientation  of  the  particles 
at  the  outflow  and  indicate  the  relative  rotation  while  traversing  the  main  body  of  the  reactor. 


Radiative  and  conductive  heating  of  the  wall  opposite  the  substrate  is  predicted  to 
be  significant.  For  a  substrate  temperature  of  773  I\  the  peak  opposite  wall  temperature 
was  predicted  to  be  541  K.  Simulations  also  predict  similar  peak  temperatures  where  the 
side  walls  join  the  opposite  wall  in  the  vicinity  of  the  substrate.  Initial  comparison  of 
these  predictions  with  thermal  imaging  data  from  the  validation  reactor  showed  reason¬ 
able  agreement  between  experiment  and  simulation.  Validation  with  other  flow  rates  and 
substrate  temperatures  is  continuing.  The  collection  of  thermal  data  is  a  collaborative 
effort  with  some  of  the  modeling  team  participating  in  the  lab  and  this  active  partici¬ 
pation  was  important  for  the  interpretation  of  the  data  and  comparison  to  simulation. 
Simulation  in  the  exhaust  portion  of  the  reactor  did  not  suggest  any  problems. 
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Figure  15:  Simulation  of  the  exhaust  portion  of  the  horizontal  reactor.  Dots  represent  traces  of  massless 
particles  as  they  traverse  this  portion  of  the  reactor. 


6  SUMMARY  REMARKS 


As  detailed  in  our  discussions  above,  early  design  of  the  HPOMCVD  reactor  dramatically 
evolved  long  before  any  physical  reactor  was  built.  The  effects  described  took  place  over 
a  period  of  approximately  2  years  and  involved  close  daily  collaboration  between  the 
material  scientists,  physicists,  and  applied  mathematicians.  Mathematical  modeling  and 
computer  simulation  were  the  essential  tools  in  the  design  effort.  The  rapid  progress  and 
substantial  savings  to  date  constitute  a  dramatic  and  powerful  argument  for  the  design 
approach  in  this  and  other  areas  of  emerging  technologies. 


ACKNOWLEDGMENTS 

This  work  was  supported  in  part  by  DOD/AFSOR  MURI  Grant  No.  F49620-95-1- 
0447  and  NASA  collaborative  agreement  NCC8-95. 


20 


References 


[1]  K.J.  Bachmann  and  G.M.  Kepler.  Heteroepitaxial  processes  at  low  and  elevated 
pressures.  SPIE  3123,  (1997)  64. 

[2]  K.J.  Bachmann,  C.  Hopfner,  N.  Sukidi,  A.E.  Miller,  C. Harris,  D.E.  A  spues.  N.A. 
Dietz,  H.T.  Tran,  S.  Beeler,  Iv.  Ito,  H.T.  Banks,  and  U.  Rossow.  Molecular  layer 
epitaxy  by  real-time  optical  process  monitoring.  Appl.  Surf.  Sci.  112,  (1997)  38-47. 

[3]  K.J.  Bachmann,  N.  Sukidi,  C.  Hopfner,  C. Harris,  N.A.  Dietz,  H.T.  Tran,  S.  Beeler, 
K.  Ito,  and  H.T.  Banks.  Real-time  monitoring  of  steady-state  pulsed  chemical  beam 
epitaxy  by  p-polarized  reflectance.  J.  Crystal  Growth  183,  (1998)  323-337. 

[4]  Y.  Bayazitoglu  and  M.N.  Ozisik,  Elements  of  Mass  Heat  Transfer ,  (McGraw  Hill, 
NY,  1988). 

[5]  T.  Bergunde,  M.  Dauelsberg,  L.  Kadinski,  Yu.N.  Makarov,  M.  Weyers,  D.  Schmitz, 
G.  Strauch,  and  H.  Jiirgensen.  Heat  transfer  and  mass  transport  in  a  multiwafer 
MOVPE  reactor:  modelling  and  experimental  studies.  J.  Crystal  Growth  170,  (1997) 
66. 

[6]  C.R.  Biber,  C.A.  Wang,  and  S.  Motakef.  Flow  regime  map  and  deposition  rate 
uniformity  in  vertical  rotating-disk  OMVPE  reactors.  J.  Crystal  Growth  123,  (1992) 
545. 

[71  R.B.  Bird,  W.E.  Stewart,  and  E.N.  Lightfoot,  Transport  Phenomena ,  (John  Wiley 
and  Sons,  NY,  1960). 

[8]  R.D.  Blevins,  Applied  Fluid  Dynamics  Handbook ,  (Krieger,  Malabar  FLA,  1984) 

[9]  D.I.  Fotiadis  and  S.  Kieda.  Transport  phenomena  in  vertical  reactors  for  metalor- 
ganic  vapor  phase  epitaxy.  J.  Crystal  Growth  102,  (1990)  441. 

[10]  A. I.  Gurary,  G.S.  Tompa,  A.G.  Thompson,  R.A.  Stall,  P.A.  Zawadzki,  and  N.E. 
Schumaker.  Mechanical,  thermal  and  flow  dynamics  issues  in  compound  semicon¬ 
ductor  MOCVD  reactor  design.  J.  Crystal  Growth  145,  (1994)  642. 

[11]  J.  Hilsenrath,  C.W.  Beckett,  W.S.  Benedict,  L.  Fanno,  H.J.  Hoge,  J.F.  Masi,  R.L. 
Nuttall,  and  Y.S.  Touloukian,  Tables  of  Thermodynamic  and  Transport  Properties , 
(Pergamon  Press,  NY',  1960). 

[12]  F.P.  Incropera  and  D.P.  DeWitt,  Introduction  to  Heat  Transfer ,  (Wiley  and  Sons, 
NY,  1985). 


21 


[13]  G.M.  Kepler,  C.  Hopfner,  J.S.  Scroggs,  and  K.J.  Bachmann.  Feasibility  of  a  vertical 
reactor  for  high  pressure  MOCVD.  Mater.  Res.  and  Eng.  B,  (submitted). 

[14]  C.A.  Larsen,  N.l.  Buchan,  S.H.  Li,  and  G.B.  Stringfellow.  Decomposition  mecha¬ 
nisms  of  trimethylgallium.  J.  Crystal  Growth  102,  (1990)  103-116. 

[15]  D.R.  Lide  and  H.V.  Kehiaian,  CRC  Handbook  of  Thermophysical  and  Thermochem¬ 
ical  Data ,  (CRC  Press,  Boca  Raton,  1994). 

[16]  R.A.  Svehla,  NASA  Technical  Report  R-132,  1962. 

[17]  D.W.  Weyburne  and  B.S.  Ahern.  Design  and  operating  considerations  for  a  water- 
cooled  close-spaced  reactant  injector  in  a  production  scale  MOCVD  reactor.  .J.  Crys¬ 
tal  Growth  170,  (1997)  77. 


22 


